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ABSTRACT 

In this paper a series of high-resolution iV-body simulations is presented in which the 
equations of motion have been changed to account for MOdified Newtonian Dynamics 
(MOND). It is shown that a low-O MONDian model with an appropriate choice for 
the normalisation as can lead to similar clustering properties at redshift z — as 
the commonly accepted (standard) ACDM model. However, such a model shows no 
significant structures at high redshift with only very few objects present beyond z > 3 
that can be readily ascribed to the low flo value adopted. The agreement with ACDM 
at redshift z = is driven by the more rapid structure evolution in MOND. Moreover, 
galaxy formation appears to be more strongly biased in MONDian cosmologies. Within 
the current implementation of MOND density profiles of gravitationally bound objects 
at z — can still be fitted by the universal NFW profile but MOND halos are less 
clumpy. 

Key words: galaxy: formation - methods: iV-body simulations - cosmology: theory 
- dark matter - large scale structure of Universe 



1 INTRODUCTION 



Although the currently favoured ACDM model has proven 
to be remarkably successful on large scales (cf. WMAP re- 
sults, Spergel et al. 2003), recent high-resolution A-body 
simulations seem to be in contradiction with observation on 
sub-galactic scales: the Cold Dark Matter "crisis" on small 
scales is far from being over. The problem with the steep 
central densities of galactic halos, for instance, is still un- 
solved as the highest resolution simulations favor a cusp 
with a logarithmic inner slope for the density profile of ap- 
proximately -1.2 (Power et al. 2003), whereas high resolu- 
tion observations of low surface brightness galaxies are best 
fit by halos with a core of constant density (de Block & 
Bosma 2002; Swaters et al. 2003). Suggested solutions to 
this include the introduction of self-interactions into colli- 
sionless A-body simulations (e.g. Spergel & Steinhardt 2000; 
Bento et al. 2000), replacing cold dark matter with warm 
dark matter (e.g. Knebe et al. 2002; Bode, Ostriker & Turok 
2001; Colin et al. 2000) or non-standard modifications to an 
otherwise unperturbed CDM power spectrum (e.g. bumpy 
power spectra; Little, Knebe & Islam 2003; tilted CDM; Bul- 
lock 2001c). Some of the problems, as for instance the over- 
abundance of satellites, can be resolved with such modifica- 
tions but none of the proposed solutions have been able to 
rectify all shortcomings of ACDM simultaneously. 

Therefore, there might be alternative solutions worthy 
of exploration, one of which is to abandon dark matter com- 
pletely and to adopt the equations of MOdified Newtonian 
Dynamics (MOND; Milgrom 1983; Bekenstein & Milgrom 



1984). It has already been shown by other authors that 
this simple idea might explain many properties of galax- 
ies without the need of non-baryonic matter (e.g. Scarpa 
2003; McGaugh & de Blok 1998; Sanders 1996; Milgrom 
1994; Begeman, Broeils & Sanders et al. 1991). MOND is 
also successful in describing the dynamics of galaxy groups 
and clusters (Sanders 1999; Milgrom 1998), globular clusters 
(Scarpa 2002) and, to a limited extent, gravitational lensing 
(Mortlock & Turner 2001; Qin & Zou 1995). A recent review 
of MOND is given by Sanders & McGaugh (2002) which also 
summarizes (most of) the successes alluded to above. 

However, there has yet to come a detailed study of the 
implications of MOND in cosmological simulations of struc- 
ture and galaxy formation, which is the aim of the current 
study. Nusser (2002) already investigated modified Newto- 
nian dynamics of the large-scale structure using the A-body 
approach. His simulations, however, are lower resolution, 
both in terms of spatial and mass resolution, making a study 
of individual objects difficult. Moreover, his implementation 
of the MOND equations is slightly different to our treat- 
ment. 

The outline of the paper is as follows. In Section 2 we 
present the way we modified our A-body code MLAPM to 
account for MOND. Section 3 introduces the cosmological 
models under investigation whereas Section 4 summarizes 
the numerical details. The analysis can be found in Sec- 
tion 5. We finish with a summary and our conclusions in 
Section 6. 
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2 THE MONDIAN EQUATIONS OF MOTION 

In an Af-body code one integrates the (comoving) equations 
of motion 



1/2 



P 



pec 

a 



P = 

which are completed by Poisson's equation 
V x • F pec (x) = -A x $(x) = -4irG(p(x) - p) 



(1) 



(2) 



In these equations x = f/a is the comoving coordi- 
nate, p the canonical momentum, V^- the divergence op- 
erator (A x the Nabla operator) with respect to x and 
Fpec(x) = — V$(a?) the peculiar acceleration field in comov- 
ing coordinates. We now need to to modify these (comoving) 
equations to account for MOND. 

Despite MOND being a modification to Newton's sec- 
ond law rather than to gravity, one has the option to actually 
interpret MOND as an alteration of the law of gravity (cf. 
Sanders & McGaugh 2002). In that case Poisson's equation 



V r • g = -4nGp(r) 
is replaced by 

V r • (pdffMl/so) <?m) = -4nGp(r) , 



(3) 



(4) 



where r = ax is the proper coordinate, go the fundamental 
acceleration of MOND and gu the MONDian acceleration 
field. Note that Eq. (3) and Eq. (4) are given in proper 
coordinates in contrast to Eq. (2) where the solution F pcc 
describes only the peculiar acceleration. 

If we now compare Eq. (3) and Eq. (4) we find the 
relation between MONDian acceleration gu and Newtonian 
acceleration g to be: 



9 = K9M/go)gM + V x h 



(5) 



The field equation (4) is non-linear and difficult to solve in 
general. But it has been shown by Bekenstein & Milgrom 
(1984) that V x h decreases like 0(r~ 3 ) with increasing 
scale. Moreover, in cases of spherical, planar, or cylindri- 
cal symmetry the curl-term vanishes. Under the assumption 
that objects forming in the Universe show at least one of 
these symmetries we are able to neglect the curl-term com- 
pletely. One might argue that this kind of symmetry is prob- 
ably very weak at high redshifts. However, later in this Sec- 
tion we are also making the assumption that MOND only 
affects peculiar accelerations (in proper coordinates) which 
are well above the MOND acceleration at early times (cf. 
Eq. (12)). We therefore presume that V x h is unimportant 
for the growth of large-scale structures as well as for the in- 
ternal properties of virialized objects (under the assumption 
that they are at least symmetric along one of their axes) . 

Using Milgrom's (1983) suggested interpolation func- 
tion 



p(x) = x(l + x 2 ) 1 ^ 2 



(6) 



one now only needs to solve 

9m - g^/gf+glf = , (7) 

to get g M as a function of g. The relevant solution of Eq. (7) 



(8) 



Eq. (8) allows us to obtain the MONDian acceleration 
gM for a given Newtonian acceleration g where gM and g 
are assumed to be parallel. 

However, we are actually solving Eq. (2) in our iV-body 
code MLAPM which gives F pcc , the Newtonian peculiar accel- 
eration in comoving coordinates. Therefore, we also need to 
derive a relation between the proper acceleration g = r, the 
proper peculiar acceleration g pec , and the peculiar acceler- 
ation in comoving coordinates F poc . The second derivative 
with respect to time of r = ax gives 



r — ax + 2ax + 'ax , 

whereas combing Eqs. (1) leads to 



(9) 



(10) 



Using the second Friedmann equation we can then rewrite 
Eq. (9) as follows 



(11) 



This equation shows that the peculiar acceleration in proper 
coordinates g pec should be defined as 



g P ec — F pcc /a 



(12) 



where F pec is the solution of Poisson's equation Eq. (2) as 
obtained by MLAPM. 

If we now make the assumption that MOND only affects 
the peculiar acceleration in proper coordinates but leaves the 
Hubble acceleration unchanged, the recipe for getting the 
MONDian peculiar accelerations in comoving coordinates 
Fm, pec is as follows: 

i) solve Eq. (2) using MLAPM which gives F pec 

ii) use Eq. (12) to transfer F pcc to g pcc — F pcc /a 2 

iii) use Eq. (8) to calculate (?m, P cc from g pec 

iv) transfer gM, pcc back to F M , P cc = a 2 g M , P ec 

v) use Fm, P cc rather than F pec in Eq. (1) 

This scheme has been employed for the simulations car- 
ried out with MLAPM described below*. We like to point out 
that our treatment of MOND agrees with the one advocated 
by Sanders (2001) in the limit (3 = 0. 



3 THE COSMOLOGICAL MODELS 

The reason for introducing MOND by Milgrom (1983) was to 
explain the fiat rotation curves of galaxies without the need 
for dark matter. Having this in mind we decided to use an 
flo value for the MOND simulation that is close to the upper 
bound allowed by Big-Bang-Nucleosynthesis and agrees with 
the recent measurements of cosmological parameters by the 
WMAP experiment (Spergel et al. 2003). Our database of 
simulations is made up of the following three runs 

* The latest release version of MLAPM includes the MOND imple- 
mentation which can be activated using -DM0ND on compile time. 
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Table 1. Model parameters. In all cases a value for the Hubble 
parameter of h = 0.7 was employed. 



label 


n 


n b 


Ac, 




CT no rm 


go [cm/s 2 ] 


ACDM 


0.30 


0.04 


0.7 


0.88 


0.88 




OCBM 


0.04 


0.04 


0.0 


0.88 


0.88 




OCBMond 


0.04 


0.04 


0.0 


0.92 


0.40 


1.2 xlO 8 



• a standard ACDM model, 

• an open, low-fio model with the same as as ACDM, 

• an open, low-fi model with MOND and adjusted a 8 . 

These runs are labeled ACDM, OCBM and OCBMond, 
respectively, and their cosmological parameters are summa- 
rized in Table 1. The OCBM model is only to be understood 
as a gauge for the MOND model run rather than an alter- 
native to ACDM. 

We can see from Eq. (12) that the peculiar accelerations 
(which are subject to MOND) are large at early times and 
therefore a "Newtonian treatment" in the early universe is 
justified. Therefore, the input power spectra to our initial 
conditions generator were calculated using the CMBFAST 
code (Seljak & Zaldarriaga 1996) with fi = fit, = 0.04 
for OCBM and OCBMond, respectively. This explains the 
choice for using the expression CBM rather than CDM. Such 
a relatively high fit value (compared to fio) actually intro- 
duces "baryon wiggles" into the primordial power spectrum 
(oscillations frozen into the plasma at the epoch of recom- 
bination which are suppressed in dark matter dominated 
models). However, fluctuations on scales of the box size em- 
ployed (B = 32/i _1 Mpc, see below) and smaller are not 
affected by it other than an overall "damping" (e.g. Eisen- 
stein & Hu 1998; Silk 1968). This damping, however, is com- 
pensated by the choice of normalisation crg orm of the power 
spectrum. We, moreover, differentiate between <7g orm and 
the actual measure of <7g =0 in the simulations at redshift 
z = because the OCBMond model requires a lower P(k)- 
normalisation <7g orm to arrive at a comparable <7g =0 value. 
This is due to a much faster growth of structures when using 
MOND which will be emphasized in more detail in Section 5. 



4 THE A/-BODY SIMULATIONS 

Using the input power spectra derived with the CMBFAST 
code we displace 128 3 particles from their initial positions 
on a regular lattice using the Zel'dovich approximation (Ef- 
stathiou, Frenk & White 1985). The box size was chosen 
to be 32k" 1 Mpc on a side. This choice guarantees proper 
treatment of the fundamental mode which will still be in 
the linear regime at z = (cf. the scale turning non-linear 
at z — is roughly 20/i _1 Mpc for the models under inves- 
tigation). The particles were evolved from redshift z = 50 
until z — with the open source adaptive mesh refinement 
code MLAPM (Knebe, Green & Binney 2001). We employed 
500 steps on the domain grid built of 256 3 cells, and in 
all three runs a force resolution of llh^ 1 kpc was reached 
in the highest density regions. The mass resolution of the 
runs is m p = 1.30 • 10 9 /i _1 M for the ACDM model and 
m p = 0.17 • 10 9 /i _1 M© for the two low-fio models, re- 
spectively. We output the particle positions and velocities 
at redshifts z = 5,3,1,0.5, and 0. These "snapshots" are 



then analysed with respect to the large-scale clustering as 
well as properties of individual objects. And even though in 
OCBM(ond) we made the assumption of the non-existence 
of dark matter we still refer to these objects as "halos"; due 
to the absence of dissipation they show similar shapes and 
density distributions when being compared to their ACDM 
counterparts as can be seen in the analysis below. 

Gravitationally bound objects are identified using the 
Bound-Density-Maxima method (BDM, Klypin & Holtzman 
1997). The BDM code identifies local overdensity peaks by 
smoothing the density field on a particular scale of the or- 
der of the force resolution. These peaks are prospective halo 
centres. For each of these halo centers we step out in (log- 
arithmically spaced) radial bins until the density reaches 
phaio(fvir) = A v i r pf, where pt is the background density. This 
defines the outer radius r v j r of the halo. However, one needs 
to carefully choose the correct virial overdensity A v i r which 
is much higher for the OCBM and OCBMond model due to 
the low fio value. The parameters used are A v i r = 340 for 
ACDM and A vir = 2200 for OCBM/OCBMond (see Gross 
1997, Appendix C, and references therein). 

Brada & Milgrom (1999) raised the issue that care 
should be taken when numerically integrating the equations 
of motion using MOND. Even more so in our case where we 
made tha assumption that the curl-term in Eq. (5) vanishes 
which gives rise to forces that are not strictly conservative. 
We therefore performed a simple test on the final output to 
assure that our implementation of MOND does not violate 
momentum conservation: the net momentum of all particles 
in the box as well as the net force should vanish. We calcu- 
lated the following sums 



E^ 



CiJ2\vi\ 



and 



^E 



F, 



(13) 



and derived values for Ci and C2 in all three models. In 
ACDM and OCBM those constants Ci,2 are of the order 
10" 4 and even though in OCBMond they are about an order 
of magnitude larger we still believe that they are sufficiently 
small. It has been shown elsewhere (Knebe, Green & Binney 
2001) that adaptive softening in general does not guarantee 
precise momentum conservation and our values for Ci and 
C2 are as expected. 



5 ANALYSIS 

5.1 Large-Scale Clustering Properties 

We start with inspecting the large-scale density field in all 
three runs. Fig. 1 shows a projection of the whole simulation 
with each individual particle grey-scaled according to the lo- 
cal density at redshift z = and z = 5. This figure indicates 
that the MOND simulation looks fairly similar to the other 
two models in terms of the locations of high density peaks 
(dark areas), filaments and the large-scale structure, respec- 
tively. One should bear in mind though that the OCBMond 
simulation was started with a much lower <7g orm normali- 
sation than the other two runs. This is in fact reflected in 
the right panel showing the density field at redshift z = 5; 
the MOND simulation is less evolved. Moreover, without 
any further analysis one might even be inclined to conclude 
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Figure 1. Comparison of the large-scale density field of the three models under investigation at rcdshift z = (left) and 2 = 5 (right). 
The upper panel shows the ACDM simulation, the middle panel the OCBMond model and the lower one OCBM. 
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1 10 100 

k [h Mpc- 1 ] 

Figure 2. Matter power spectra for all three model at rcdshift 
z = (thick lines) and 2 = 5 (thin lines). 

from Fig. 1 that the MOND model is more strongly clus- 
tered at z = which is affirmed by the slightly higher erf 
value given in Table f. But as we will see later on this is 
not necessarily true; we can confirm a higher amplitude of 
the two-point correlation function for galaxies in OCBMond 
(cf. Fig. 5) while at the same time showing a comparable 
amplitude in the matter power spectrum on small scales. 
The latter can be viewed in Fig. 2 where we plot the matter 
power spectrum for all three models at redshift z — 0. There 
we observe that the OCBMond model shows a slightly larger 
amplitude for k < 0.8hMpc -1 (scales close to the fundamen- 
tal mode). Nusser (2002), who treated the MOND equation 

similarly (but differently) to us^, already pointed out that 
the linear evolution of the growing mode solution for the 
density contrast S scales like 8 oc a 2 as opposed to New- 
tonian theory where S oc a. This explains why the OCB- 
Mond model with the (initially) low <7g orm normalisation 
outruns the evolution of the Newtonian OCBM simulation. 
In other words, the MOND model had to be started with 
a lower as normalisation to provide competitive results at 
redshift z = 0. Sanders (2001) also noted that due to a 
much faster growth of structures in MOND universes the 
amplitude of P(k) for purely baryonic models matches the 
standard ACDM cosmology. We also like to point out that 
the required value <Tg orm = 0.4 needed to bring the MOND 
simulation into agreement with the standard ACDM model 
is closer to the COBE normalisation <jg OBE £s 0.1 for that 
particular cosmology. Another feature worth noting is that 
the OCBMond power spectrum does not show a distinctive 
"break" due to the transfer of power from large to small 
scales as seen in both the ACDM and the OCBM model. 

t Nusser (2002) did not use Milgrom's interpolation function as 
given by Eq. (6) but rather a spontaneous transition from New- 
tonian to MOND accelerations. 




10 10 10" 10 12 10 13 10 14 10 10 10" 10 12 10 13 10 14 



Figure 3. Cumulative mass function for all three models at red- 
shifts z = 5, 3, 1 and 0. 

In Fig. 3 we are plotting the cumulative mass function 
of objects identified using the BDM technique (Klypin & 
Holtzmann 1997). This figure highlights again that the hier- 
archical formation of gravitationally bound objects is driven 
much faster in the MOND simulation but being initiated 
at later times. At a redshift of z — 5 we can see that the 
abundance of objects on all mass scales is nearly identical 
in the OCBM and ACDM model with a much lower am- 
plitude for OCBMond. Whereas the evolution for the New- 
tonian OCBM run already ceases at a redshift of around 
z s top — 1 /Ho — I ~ 24 we still see a very strong increase (by 
more than one order of magnitude) in the number density of 
objects of all masses in the OCBMond simulation. To em- 
phasize on this, Fig. 4 shows the (integral) abundance evo- 
lution of objects with mass M > W 11 ^ 1 M Q . The OCBM 
model undoubly experiences very little evolution from z ~ 5 
to today whereas both other models show a very steep evolu- 
tion. The discrepancy between the OCBMond and the other 
two models at redshift z = 5 can, however, be ascribed to 
the lower initial <7g orm value though; OCBMond was set up 
with much smaller initial density perturbation which only 
grew to a comparable level of clustering via the effects of 
MOND. This is again in agreement with the findings of 
Sanders (2001) who showed that the collapse of spherically 
symmetric overdensities becomes MOND dominated for red- 
shifts z < 5 and hence starts to outrun Newtonian models 
(cf. Fig.fTin Sanders 2001). 

The question now arises to what degree the (formation) 
sites of halos in ACDM and the two OCBM models are corre- 
lated. To this extent we calculated the two-point correlation 
function of the 500 most massive objects and the result can 
be found in Fig. 5. We chose to use a fixed number of halos 
rather than a mass cut for the calculation of £ ga i (r) in order 
not to introduce an artificial bias. We do have far more ob- 
jects of a given mass in the ACDM model and therefore a 
mass limit M m i n used with the estimator for £haio(f) would 
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Figure 5. Two-point correlation function at rcdshift z = for 
the 500 most massive objects. 



lead to different correlation amplitudes. The agreement be- 
tween the Newtonian OCBM and the ACDM model is not 
surprising. As already pointed out by other authors, the cor- 
relation function is expected to be (nearly) identical in cases 
of equal as normalizations, irrespective of the cosmological 
model (Martel & Matzner 2000). Moreover, if the model is 
fixed and only the as normalization varied it should leave no 
imprint on £ ga i(r) (Croft & Efstathiou 1994). But the OCB- 
Mond model stands out again having a much higher ampli- 
tude on all scales. This clearly attributes for the differences 
already mentioned in the discussion of Fig. 1 and indicates 
that the OCBMond model is more evolved at redhsift z — 
than the other two models even though the matter power 
spectrum has a lower amplitude on corresponding scales: 
structure formation in MONDian cosmologies is even more 
biased than in Newtonian models. This is in agreement with 
findings that small scales enter the MOND regime before 
large-scale fluctuations (cf. Nusser 2002; Sanders 2001). 



Figure 6. Most massive galactic halo in ACDM (upper panel), 
OCBMond (middle panel) and OCBM (lower panel). The line 
in the lower right corner of each panel is a scale indicating 
100ft- 1 kpc. 



5.2 Galactic Halos 

Having analysed the large-scale clustering properties we now 
turn to the investigation of the internal properties of galactic 
halos. To this extent we use the halo catalogues based on 
the BDM code (Klypin & Holtzmann 1997) and neglecting 
objects less massive than M < W 11 ^ 1 M Q again. This sets 
the minimum number of particles per halo to 77 for ACDM 
and 488 for OCBM(ond). 

A visualization of the density fields throughout the most 
massive BDM halo is given in Fig. 6. It is quite striking that 
neither of the low-f2o halos shows substantial substructure. 
However, this is easily understood for OCBM, because in 
a low-density universe structure formation ceases at early 
times (zstop — l/^o — !)• This means that clusters in such 
cosmologies should show fewer substructure since they had 
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more time to virialize (cf. Knebe & Miiller 2000). But this 
explanation obviously does not hold for OCBMond as nearly 
all halos in that particular model formed exceptionally late 
(cf. Fig. 4). 

The most interesting question, however, is probably the 
shape of the density profile and the rotation curve for halos 
in MONDian cosmologies, respectively. Fig. 7 now shows the 
matter profile of the most massive halo in all three models 
along with fits (thin solid lines) to a Navarro, Frenk & White 
(NFW) profile (Navarro, Frenk & White 1997) 

PnfwWoc ^^I^^ . (14) 

The scale radius r s is being used to define the concentration 
of the halo c = r v i T /r s where r v i r is the radius at which 
the density reaches the virial overdensity A v i r w 340 and 
A vir « 2200 for ACDM and OCBM(ond), respectively. 

We observe that even for the OCBMond model the data 
is equally well described by the functional form of a NFW 
profile (out to the virial radius, at least). However, the cen- 
tral density of that halo in the OCBMond model is lower 
than in ACDM and especially in OCBM. The high central 
density for OCBM is readily explained by the fact that ha- 
los in that particular cosmology form at a time when the 
universe is still very dense. This result is also supported by 
the values of the concentration parameter presented in Ta- 
ble 2: the most massive object in the MOND model shows 
the lowest concentration c, mostly due to the late onset of 
formation as observed in Fig. 4. When inverting the density 
profile into a (Newtonian) circular velocity curve by simply 
using «circ( r ) = GM(< r)/r, this also entails a shift of the 
maximum of the rotation curve to higher radii as can be seen 



in Fig. 8 



2r s , Navarro, Frenk & White 1997). The 



functional shape of the rotation curves for an NFW density 
profile is given by 



^circ 
o,2 



1 ln(l + cx) - (cx)/(l + cx) 
x In (1 + c) - c/(l + c) 



(15) 



with w v i r being the (Newtonian) circular velocity at the virial 
radius r vir and x — r/r vil . The drop of the maximum of the 
rotation curve by about a factor of 3.2 is merely a reflection 
of the scatter in mass for the most massive halo throughout 
the three models. As can be seen in Table 2 the halo is more 
than ten times as massive in ACDM than in OCBMond. This 
should give an about 2.7 times higher «mix as the scaling 
between those two quantities is roughly f^lS oc M v Y r 3 . This 
scaling relation can be derived when using i™ ~ 2/c (cf. 
Bullock et al. 2001b and Navarro, Frenk & White 1997) 
together with Eq. (15) giving 



ln(l + c)-c/(l+c) 



(16) 



If we furthermore assume c oc M -013 as shown by Bul- 
lock et al. (2001b, cf. Eq(18)) we find that the ^TT-factor in 
Eq. (16) is roughly constant for the mass range under con- 
sideration. And as v v i T oc M 1 ^ 3 (which simply follows from 
v 2 oc M/r and r oc M 1//3 for a spherical overdensity) the 
same scaling then holds for u™ r a c x explaining the observed 
drop of the maximum of the rotation curve for the OCB- 
Mond model. 

However, for the OCBMond model we should take into 
account that accelerations are not Newtonian and hence 
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Figure 7. Density profile for the most massive halo in all three 
models. The vertical thin lines are indicating the virial radii. 
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Figure 8. Rotation curve of the most massive halo. 



f<Lc( r ) = GM(<r)/r does not hold. Therefore we also show 
in Fig. 8 the actual MONDian rotation curve calculated as 
follows: the Newtonian acceleration g — v 2 /r is transfered 
to the MONDian acceleration qm according to Eq. (8). gtA 
in turn is used to calculate vm(t) = y/gur. The resulting 
Wcirc(r) is labelled "MOND" in Fig. 8. We note that the 
MONDian velocities are actually larger than the Newtonian 
ones bringing them closer to the ACDM model. This has 
implications for dynamical mass estimates of galaxy clus- 
ters as presented in Sanders (1999). Sanders showed that 
the dynamically estimated mass of galaxy clusters is too 
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Figure 9. Acceleration curve of the most massive halo. The solid 
vertical line indicates the MONDian accleration parameter go = 
1.2 X 10 8 cm/sec 2 . 



large compared to the observed mass when using Newto- 
nian physics. However, MOND rectifies dynamical masses 
bringing them into better agreement with observed masses. 
And a similar phenomenon can be observed in Fig. 8: the 
MONDian curve for OCBMond would be measured obser- 
vationally and hence be translated into a too high cluster 
mass of M vir w 2.8- lO 13 ^ 1 M when using Newtonian dy- 
namics. MOND, on the contrary, would give the real value 
of M vir = 0.3 • lO 13 ^ 1 M Q . Later on we will see that this 
leaves an imprint on the (radial) distribution of the spin 
parameter, too. 

In addition to the rotation curve we also show the ac- 
celeration as a function of radius in Fig. 9. The object that 
formed in OCBMond has MONDian accelerations through- 
out the whole halo whereas the central parts of the objects 
in the other two models are above the universal acceleration 
go indicated by the thin solid line. 

We like to stress that in both Figs 7 and 8 the New- 
tonian curves for the OCBMond model are only plotted 
for completeness; they do not carry observable informa- 
tion as the halo actually follows MONDian physics rather 
than Newtonian. However, they are educational in a sense 
to gauge the importance MOND has on the internal struc- 
ture of the halo. 

Table 2 now lists some internal properties in addition to 
the ones already mentioned, i.e. the velocity dispersion a v , 
the virial radius r v i r , the triaxiality T, ellipticities ci and e2, 
the spin parameter A as well as best-fit parameters Ao and 
a\ when fitting the probability distribution P(X) to a log- 
normal distribution. The spin parameter A was calculated 
using the definition given in Bullock et al. (2001a) 

J 



A = — , (17) 

V2M vir «virrvir 

which proved to be a more stable measurement than the 
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Figure 10. Spin parameter distribution in all three models. Lines 
show fits obtained using the log-normal distribution given by 
Eq. (18). 



usual A' = Jyf\E\/(GM 5/2 ) definition. For the OCBMond 
model w v i r is again the MONDian circular velocity at the 
virial radius. The distribution of A, P(A), has been fit- 
ted to a log-normal distribution (e.g. Frenk et al. 1988; 
Warren et al. 1992; Cole & Lacey 1996; Mailer, Dekel & 
Somerville 2002; Gardner 2001) 



P(A) 



\V2tt ao 



exp 



ln 2 (A/A ) 
2a 2 



(18) 



The results are presented in Fig. 10 showing that the OCB- 
Mond distribution P(A) is nearly indistinguishable from the 
other two models. However, the reduced \ 2 value for OCB- 
Mond is about a factor of two larger (cf. Table 2) . As been 
noted by Mailer, Dekel & Sommerville (2001) the log-normal 
distribution given by Eq. (18) is not as good a fit to models 
with halos that are recent merger remnants. This is defi- 
nitely one of the effects that has an influence on the spin 
parameter distribution for the OCBMond model as we ex- 
pect a high level of recent merger activity (cf. Fig. 4). 

We also calculated the radial distribution of A(< r) 
throughout the halos and present the result for the most 
massive one in Fig. 11. We see that A(< r) is roughly 
constant for the Newtonian models of structure formation 
whereas there is a sharp increase of A(< r) in the MOND 
halo towards the virial radius. This implies that the mate- 
rial in that particular halo moves on more circular orbits and 
the object itself is closer to solid body rotation, respectively. 
This result is in agreement with our previous finding where 
we showed that the MONDian rotation curve is not given by 
simply "inverting the density profile" as in the Newtonian 
case (cf. Fig. 9); the velocity in the outskirts of the halo is 
much larger and hence we expect a rise in A(r), too. 

We are now going to focus on the shape of the halos. 
Firstly, we show measurements of the overall shape of halos 
and to this extent we calculated the eigenvalues a > b > c of 
the inertia tensor. They were in turn used to construct the 
triaxiality parameter (e.g. Franx, Illingworth & de Zeeuw 
1991): 



© 0000 RAS, MNRAS 000, 000-000 



Galactic Halos in MONDian Cosmological Simulations 9 



0.05 



0.04 



0.03 



0.02 



. . .ULbMOnO ■-■ 










: OCBM r .J-" 




L 






r ACDM r' 


















0.10 










b 


, , 1 








0.01 


100 


1000 



r [h ? kpc] 

Figure 11. Radial distribution of spin parameter A. The vertical 
lines are again showing the virial radius. 



Table 3. Mean triaxiality parameter (T) and mean ellipticities 
(ei)> ( e 2) when averaging over halos with M > lO 11 ^. -1 Mg. 



label 



(T) (ei) <e 2 ) 



ACDM 0.55 0.25 0.13 

OCBMond 0.61 0.28 0.17 
OCBM 0.55 0.21 0.12 



T = 



(19) 



The triaxialities T are accompanied by the ellipticities 
b 



ei = 1- 



e 2 = l- 



(20) 



and the values for the most massive object are summarized 
in Table 2. The mean values (T), (ei), and (e2) when aver- 
aging over all halos more massive than W 11 ^ 1 Mq can be 
found in Table 3. We observe a trend for the MOND halos 
to be more triaxial with higher ellipticities at the same time. 

Secondly, we like to quantify subclumps within the virial 
radius of the halo itself. One possibility to measure the sub- 
structure content of a halo is to calculate the radial profile 
of the density dispersion 



1 



N(r) 



N(r) 

£ 

i=l 



Mr) - (p(r)) 
(P(r)) 



(21) 



where pi(r) is the local density at a particle position which 
was estimated using the nearest 20 neighbors and (p(r)) the 
average taken over all N(r) particles within a spherical shell 
[r, r + dr] . We used the same binning as already applied 
to the density profile and the rotation curve, respectively. 
The result for the most massive halo in all three runs is 
presented in Fig. 12. This figure shows that the dispersion 
is always smallest for the MOND halo, at least out to the 
virial radius: the density at each individual particle position 
is always close to the mean density. If there were subclumps 
present one would expect to find peaks in the af (r) curve 
due to local deviations from the mean density profile. 
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Figure 12. Variance <r^ (r) for the most massive halo. The vertical 
lines are indicating the virial radius of the respective halo. 



6 SUMMARY AND DISCUSSION 

In this paper we presented three cosmological simulations: a 
fiducial standard ACDM model, a very low-no model, and 
the same Qo model but with MONDian equations of mo- 
tions. Putting aside the classical arguments against MOND, 
we showed that structures found in a MONDian low-f2 uni- 
verse are not significantly different from the standard ACDM 
model. However, we derived several differences which might 
easily validate or rule out MONDian cosmologies observa- 
tionally. Our main results can be summarized as follows: 

• even though the OCBMond run was set up using a 
low value for <rg orm = 0.4 the simulation finished with 
(j| =0 = 0.92 (which is close to the cluster normalisation 
of the ACDM model), 

• the OCBMond model shows an extremely fast evolution 
of the number density of halos for redshifts z < 5, 

• virialized objects in a MONDian universe are slightly 
more correlated, 

• the density profile of MOND objects still follows the 
universal NFW density profile, but they are 

• far less concentrated, 

• show less substructure, and 

• are closer to solid body rotation. 

We conclude that the most distinctive feature of a MON- 
Dian universe is the late epoch of galaxy formation; in a 
MONDian universe that is normalised to match a ACDM 
model at redshift z — we expect to not observe any galax- 
ies until recently (z < 5). This actually holds for any low-as 
universe either MONDian or Newtonian, but only the as- 
sumption of MOND can bring such a model into agreement 
with observations of the local universe again. Another in- 
teresting finding is that the outer parts of MOND halos are 
closer to solid body rotation than their standard ACDM 
counterparts even though the overall distribution P(X) is 
nearly indistinguishable from the ACDM model. 

However, there have still been many assumptions made 
during the course of this study which are hard to justify 
and hence all results are to be understood simply as prelim- 
inary until our understanding of MOND actually improves. 
We not only assumed that MOND only affects peculiar ac- 
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Table 2. Properties of the most massive gravitionally bound object. The mass M is measured in h 1 Mq, velocities in km/sec and radii 
in h^ 1 kpc. x 2 is actually the reduced \ 2 value as returned by IDL's curvef it routine. 

model M a v v v i T r v i r r s c T e\ e2 A Ao u x X 2 

ACDM 5.6-10 13 657 553 790 91 8.7 0.65 0.30 0.18 0.025 0.040 0.53 0.8 

OCBMond 0.3-10 13 286 473 300 59 5.1 0.81 0.30 0.23 0.057 0.034 0.52 1.5 
OCBM 0.8-10 13 373 271 418 27 15.5 0.80 0.21 0.16 0.023 0.034 0.52 0.9 



celerations in proper coordinates but also neglected the curl- 
term in Eq. (5). The effect of this rotational component is 
not clear, but it is guaranteed to decrease rapidly on large 
scales and vanish completely for objects that have spheri- 
cal, planar or cylindrical symmetry. Another disclaimer is 
that MOND is based on the non-existence of dark matter, 
but our simulations only model gravity; if the universe con- 
sist only of baryons one definitely would need to include gas 
physics to be able to make credible predictions for inter- 
nal properties of galaxies. But this is beyond the scope of 
this study. Nonetheless, we have shown that cosmology with 
MOND does not necessarily lead to completely odd results. 
Our treament of MOND though gives clustering properties 
comparable to the favorite concordance ACDM model. 
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